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1. Introduction 



We tackle the problem of subpixel object detection in image sequences which arises for in- 
stance in infrared search and track (IRST) applications. In this context, the target signature 
is proportional to: 



/•i+0.5 /■j+0.5 

s e [i,j]= / / h (u-e 1 ,v-e 2 )dudv. (1) 

Ji-0.5 Jj-0.5 

s e [i,j] represents the percentage of light intensity at pixel e = (ei, e 2 ) refers to the 
object random subpixel position and h a is the optical point spread function (PSF). Ac- 
cording to common sensor design, the energy of the signal component s = as e is almost 
concentrated on a single pixel. However, contrary to the amplitude a which is unknown too, 
dependence on the location parameter e is highly nonlinear. Its influence in our application 
is rather significant because of aliasing and unless a velocity model is available, object sub- 
pixel position is hardly predictable from frame to frame. Actually, common sensor design 
leads to an image spot downsampled by almost a factor 5. We can see on Figured] the en- 
ergy loss at central pixel according to subpixel location and the random change in spatial 
pattern due to aliasing. This phenomenon has a major impact on detection performance as 
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Fig. 1 . Examples of image spots for different cross-marked subpixel positions (win- 
dows of size 5x5 pixels). Sensor design parameter r c is set to its common value of 
2.44 (see section HJ). 



shown thereafter. To our knowledge, this pitfall has not been addressed yet in the litera- 
ture. A prevailing opinion stands that there is no signature information in subpixel objects. 
Indeed, the different authors dealing with small object detection have concentrated on clut- 
ter removing, 13 multi- or hyper- spectral fusion 4 5 and multiframe tracking methods. 6-8 We 
focus here on the processing of a single frame. In section [2l we formulate the detection 
problem in the classical model of a signal in additive Gaussian noise [9, ch.2-4]. When the 
signal is deterministic, Neyman-Pearson strategy yields the conventional matched filter. In 
the present case, the signal from the target depends on unknown parameters and we have 
to deal with a composite hypothesis test. A common procedure is given by the generalized 
likelihood ratio test. But the " nuisance " parameters a and e can also be considered as ran- 
dom variables with known distributions (some a priori density functions in the Bayesian 
terminology), then the straightforward extension of the likelihood ratio test is to integrate 
the conditional distribution over a and e. When modelling the signal component as a sam- 
ple function, we could also think of the class of random signal in noise detection problems, 
which have essentially been studied in the Gaussian case. Unfortunately, considering s e as 
a random vector, its empirical distribution proves to be highly non Gaussian when e is uni- 
formly sampled. For instance, the histogram of the central pixel depicted on Figure|2]shows 
that a Gaussian fit is not satisfactory at all. In section|3l we define more precisely the optical 
system model used in our numerical experiments. We consider both a Gaussian white noise 
and a fractal noise of unknown correlations generated by a standard technique of spectral 
synthesis. Section|4]is devoted to the position estimation problem, i.e. estimation of param- 
eter e. We propose two estimators that take into account the fact that the signal amplitude 
a is also unknown. We demonstrate the performance of these estimators in terms of mean 
square errors. As for the detection problem, we finally illustrate the expected improvement 
in quality brought by a correctly sampled optics compared to common sensor design. 
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Fig. 2. Empirical distribution of the image-spot central pixel s e [0, 0] for a uniformly 
random position e ~ W[-o.5,o.5[ 2 - 



2. Detection problem 

We consider a local detection window sliding across the image. The problem is to decide 
whether an object is present or not at the window central pixel. This is a binary test which 
typically reads as follows: 

H : z = n 
Hi : z — as e + n 

where z is the vector collecting the window data, s = as e is the object response (signal 
vector) and n the additive Gaussian noise. The signature shape is known and deterministic, 
so that s only depends on the two unknown parameters a G R and e e £ = [—0.5, 0.5 [ 2 . 
The noise vector n is supposed to be centered (in practice we first remove the empirical 
mean from the data) with a known or previously estimated covariance matrix R. Thus, if we 
assume that n is independent from s, the following conditional distributions are Gaussian: 

p(z|H ) ~ Af(0,R) 
p(z|Hi,a,e) ~ Af(as e ,R) 

Let first assume that parameters a and e are given. The problem amounts to a simple 
hypothesis test which is to detect a deterministic signal in a Gaussian noise. The Neyman- 
Pearson strategy or likelihood ratio test (LRT) is given by: 



Hi 

p(z|Hi,q,e) > 

p(z|H ) < 
H n 



threshold (4) 



It is equivalent to classical matched filtering which simply compares the statistic aT e (z) = 
a s l e R~ l z with some threshold. 

A. Pixel matched filtering 

The exact object location being unknown in practice, we could assume by default that 
e = e = [0,0], i.e. the object is at the center of the pixel, whereas the true location 
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would correspond to e = e*. Thus, the detector which consists in thresholding the pixel 
matched filter (PMF) aT eo (z) is optimum provided e* = e . Otherwise it is mismatched 
and therefore suboptimum. Since conditional distributions of T eo (z) under each assump- 
tion are Gaussian, we easily get the expression of the probalility of detection Pd and of 
false alarm Pf a . Corresponding receiver operating characteristic (ROC) curves for critical 
values of e* are depicted on Figure [3] They clearly show that the PMF performances are 
significantly worse as e differs from e*. But beyond extreme situations (related to a true 
target location between two or four pixels instead of the center), the "mean curve" repre- 
sents the average statistics over uniformly random positions. Compared to the ideal curve, 
we can see that the price paid if one neglects the random location is rather high even at 
favorable signal-to-noise ratio. For a SNR of 15dB and at a Pj a of 1CT 4 , probability of 
detection decreases from nearly 1 to 0.8. 
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Fig. 3. Examples of pixel matched filter theoretical ROC curves for e* = e in blue 
(ideal curve), e* ^ e in green and red (worst case), and finally the mean curve in 
black for uniformly sampled e* (SNR = 15dB). 



The object response also depends (linearly this time) on the amplitude a, which is 
generally unknown. Yet, assuming strictly positive amplitude, we see that whatever a > 
0, thresholding aT eQ (z) gives the same ROC curve as thresholding T eo (z). Without any 
assumption on a, a classical solution is to estimate it by maximum likelihood (ML). Indeed 
under the Gaussian noise assumption, the optimum in a for a given e is explicit: 

d(e) = argmax p(z|H 1; a, e) 

= argmin {(2 — as e ) t R~ 1 (z — as £ )} 

siR^z 

and then the " generalized " pixel matched filter (referred to as GPMF) is equal to 

Is* R^zl 2 

a(e )T eo (z) = ' 60 ' . (6) 
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B. Subpixel detectors 

Our aim is to build refined detectors that improve performance of the above GPMF in taking 
into account the variability of the object signature due to its random subpixel location. 
Several solutions may be used. We first recall the most popular one. 

1. Generalized likelihood ratio test 

ML estimation of the two unknown parameters leads to the generalized likelihood ratio test 
(GLRT): 

max (ai£) p(z|Hi,q,e) 
' M) " p(*|Ho) 

(7) 

= > threshold. 

It consists in estimating the amplitude a and the possible object location e by computing: 

e ML = argmax p(z|Hi, a(e), e) 

e<=£ 



arg max 



s{R- 1 z\^ (8) 



then thresholding the estimated filter a ML % ML (z) where d ML = a(e ML ) is given by equation 
©: 

& ML % ML (z) = 1 e - . (9) 



si R-'s, 



2. Exact likelihood ratio test 



In a Bayesian approach, we propose to consider the two unknown parameters a and e as 
realizations of independent random variables with given probability density functions p(a) 
and p(e). Then the optimal procedure is the exact likelihood ratio test (ELRT). 

To compute the density function of data under Hx and to get the likelihood ratio, we 
have to integrate the conditional density p(z|Hi, a, e) over prior distributions of the nui- 
sance random parameters a and e. The likelihood ratio can be expressed as: 

= p(z|Hi) = f e / R p(z|Hi, a, e)p(at)p(e) da de 
[Z) p{z\R ) " p(z\U ) 

Given prior distributions p(a>) and p(e), C(z) is the optimal Neyman-Pearson test whenever 
a and e really satisfy the models p(a) and p(e) . By default we choose a " non-informative " 
prior for a and we adopt a uniform distribution inside the pixel for e, which seems to be 
quite a reasonable assumption for the subpixel target position. So we get: 

C(z) ex / — 1 exp \ i-g- ^ de, (11) 

Unfortunately, because of intricate nonlinear dependence of s e on e, explicit integra- 
tion over e appears to be not tractable and probability distribution of C(z) is not as simple 
as the one of T €Q (z). A quadrature approximation is required to compute C(z) whereas 
derivation of its density requires Monte-Carlo simulations. 



5 



3. Approximate likelihood ratio test 

In equation (fTTT) . the double integral over e can be approximated up to any desired accuracy 
using some quadrature rule and evaluating the integrand f(e\z) at discrete samples 6 
£ = [—0.5, 0.5 [ 2 . But, for sake of computational efficiency, we propose to use a coarsest 
approximation of the likelihood ratio (ALRT) based on a bidimensional trapezoidal rule 
which only involves the 9 half -pixel positions. 

4. Subspace model 

One alternative to this probabilistic viewpoint can be built on a geometric approach that 
restricts the signal vector s = as e to vary in some P-dimensional subspace, with P lower 
than the vector size. 10 The observed data under H4 are rewritten as: 

p 

z ~ Sa + n = a p s p + n, (12) 
p=i 

where the structural matrix S is formed by P independent vectors s p . Coefficients a p of 
the linear combination are the new parameters that describe the signal variability. Thanks 
to linearity, ML estimation of vector a has an explicit solution (which is identical to the 
least squares estimator): 

a ML = (S'R^S^S'R^z (13) 
and GLRT amounts to threshold the following statistic: 

V{z) = z t Br 1 S(S t Rr 1 S)- 1 S t Rr 1 z. (14) 

Matrix S only depends on e, a being a scale parameter. In practice, it is identified by 
discretizing £ , making a singular value decomposition and retaining the singular vectors 
s p corresponding to the P greatest singular values. We choose P = 1 which gives better 
results than higher orders. Therefore under hypothesis Hi, z ~ aiSi + n and V(z) is 
identical to GPMF with s 1 replacing s eo . 

3. Application to optical imagery 

A. Optical system 

In our application, we can model the imaging system by a diffraction-limited, unaberrated 
optics with circular aperture and incoherent illumination. 11 ' 12 The object signal pattern s e 
is then given by the integration of h a on each pixel (see equation d]), where h is the radial 
point spread function (PSF) defined by the Airy disk: 

2 

, p = Vm 2 + v 2 . (15) 

J\ is the Bessel function of the first kind and r c = v c j v s designates the normalized cut- 
off frequency (y s is the sampling frequency and v c = D/\is the radial cut-off frequency 
defined by the ratio of the lens aperture diameter D over the wavelength A). FigureHJdepicts 
the two-dimensional PSF and a slice along one diameter, as well as their Fourier transform. 



h Q (u,v) 



TV 



Ji{-Kpr c ) 
P 
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Common sensor design uses r c = 2.44 so that the pixel size is equal to the width of the 
main lobe of the PSF However, this implies a downsampling factor of v n /v s = 2r c = 
4.88 (where v n = 2v c is the Nyquist frequency). In the following section, we present 
some numerical results of detection performance considering this classical sensor design. 
Examples of image spots s e have been represented on Figured] for different values of e. 

Remark 1 We have the following property: 




1. 



P 



Fig. 4. Left: radial point spread function h (u,v) on the top and slice along a di- 
ameter on the bottom. Right: corresponding optical transfer function h (u u , u v ) and 
slice along a diameter (r c = 2.44). 



B. Numerical results 

Performances of the five classes of detectors have been compared in terms of ROC curves: 
GPMF, GLRT on a and e, ELRT, ALRT and finally GLRT with the subspace model de- 
noted SM-GLRT. Probabilities of detection and false alarm are deduced from the empirical 
distributions of these statistics under each hypothesis by generating samples of Gaussian 
noise n and uniformly distributed e in £ = [—0.5, 0.5 [ 2 . The amplitude is assumed to be 
unknown but set to a constant value a in the simulations since we have no information 
about a reliable prior distribution p(a). 

We first consider the case of a Gaussian white noise n ~ A/"(0, a 2 ) . The signal-to-noise 
ratio is then defined by: 

fa 2 E\ 

SNR=10 1og 10 (^J (16) 



with E= I (sel*,]}) 2 de 
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Fig. 5. Empirical ROC curves in the Gaussian white noise case with common sensor 
design (r c = 2.44). 



For common sensor design (r c = 2.44), the average energy of the image spot is E ~ 0.52. 
The ROC curves are depicted on Figure \5\ for two different SNR. It shows that the GLRT, 
the ELRT (actually, a refined approximation of it) and the coarse approximation ALRT 
give significantly better performance than the SM-GLRT and GPMF. We also see that the 
performance gain is greater for high SNR whereas it tends to be rather small for low SNR 
and low probability of false alarm. Conversely, if the latter detectors are computationally 
cheap, including the ALRT, this is not so for the GLRT and the ELRT, which are much 
more intensive. 

As complementary results, we have tested the five detectors on a fractal background 
image generated by a variant of the ppmforge software 0. The synthesis algorithm depends 
on the auto-similarity parameter H called Hurst parameter and which is set to 0.7 in this 
experiment. The resulting image depicted on Figure [6] is a realistic simulation of a cloud 
scene. The covariance matrix R of this stationary background is estimated by the empir- 
ical correlations on the whole image. We then compute the performance of the different 
detectors for a given target amplitude as illustrated on Figure [7J The ROC curves look quite 

'http://h30097.www3.hp.eom/demos/ossc/man-html/manl/ppmforge.l.html 
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different from the white noise case but we notice again that the GLRT, ELRT and ALRT 
have similar performance and provide a significant detection gain in comparison with the 
GPMF or the SM-GLRT. 
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Fig. 6. Simulation of a cloud fractal image of 200 x 200 pixels (Hurst parameter 
H = 0.7). 
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Fig. 7. Empirical ROC curves obtained on the fractal image of Figure [6] for a true 
(but assumed unknown) target amplitude a = 60 gray levels. 



C. Influence of the optics 

Besides the perfecting and evaluation of subpixel detectors, one additional motivation of 
this work is to analyze the influence of aliasing on detection performance. This is the reason 
why we have also tested the detectors on a correctly sampled optics in order to compare 
their performances with those obtained using a common sensor design. In the correctly 
sampled design, the focal plane is sampled at Nyquist frequency (implying a denser sensor 
array or a smaller lens diameter) so that aliasing is suppressed. Parameter r c of the PSF 
is equal to 0.5 and the signal energy is now spread over several pixels (E ~ 0.08). By 
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Fig. 8. Examples of image spots corresponding to a correctly sampled optics (r c = 
0.5) to be compared to those of Figured] 



comparison, Figure [8]presents the examples of image spots corresponding to such a design. 
Detection performances are depicted on Figure [9] on the right for a SNR of 15dB. We see 
that improved detection has just a moderate impact in this situation. The five detectors have 
a quite similar behavior but at the same SNR they perform much better than in the aliased 
case. The gain in Pf a amounts at least to a factor 10 for all the detectors. Such a result 
speaks in favour of using a denser focal plane for point target detection. 




Fig. 9. Empirical ROC curves in the Gaussian white noise case for a same SNR 
= 15dB, with common sensor design on the left (r c = 2.44) compared to a correctly 
sampled optics on the right (r c = 0.5). 



4. Performance of subpixel position estimators 

Up to now we have focused on the detection strategy. In a second step, once a potential tar- 
get is detected on a given pixel, we are also interested in accurate estimation of its subpixel 
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position. Such a problem has already been addressed, in particular for star position estima- 
tion in astronomical applications. 13 Several types of estimators are possible. We consider 
here the maximum likelihood (ML) estimator and following the Bayesian approach intro- 
duced previously the posterior mean (PM). It is important to note that the signal amplitude 
a is also unknown and therefore we have to estimate it or integrate over it. Indeed it is not 
valid to suppose that the amplitude is known in the context of IRST. 

The ML estimator of e is given in equation ([8]) by replacing a with its estimate e. Actu- 
ally, e ML and a ML = d(e ML ) are identical to joint maximum a posteriori (MAP) estimators 
with non-informative priors on the two parameters. 

The PM estimator is defined as: 



J ep(e\E 1 ,z)de (17) 



where the posterior law is deduced from Bayes'rule: 



p(e\H u z) = ri«P';«)*> (18) 
p(*|Hi) 



p( e ) f 

— — / p(;z|Hi, a, e)p(a) da. 
J R 



p(z 

So, we have to integrate over a and then over e. As previously we consider a diffuse a 
priori on R for a and a uniform law on £ for e. We get the following expression in the 
same way as for the likelihood ratio in equation ([111) : 

1 ( I s* 7? _1 zl 2 1 

p(c|Hi, z) oc - exp ' \ — j 1 . (19) 



We have studied the performance of these two estimators in terms of average mean 
square error (MSE). In practice, the optimization or the integration over e are approxi- 
mated numerically by considering a finite discrete grid of 20 x 20 values e k E 8. Given a 
true position e*, bias and variance of an estimator e are estimated thanks to Monte-Carlo 
simulations. We consider the case of a Gaussian white noise and we vary the signal-to- 
noise ratio. Figure [10] on the left compares ML and PM estimators to the pixel estimator 
which assumes by default that the target location is at the center of the pixel (e = (0, 0)) 
and whose MSE is equal to 1/12. At favorable SNR, the two subpixel estimators are far 
better than the default estimator but the gain decreases when the noise becomes important. 
For a SNR of 15dB, the ML yields an error similar to the default estimator while the PM 
notably has a twice smaller error. By comparison, Figure[lO]on the right shows the estima- 
tion performances obtained in the unaliased case (r c = 0.5) for equivalent signal-to-noise 
ratios. ML and PM logically perform better since the signal is correctly sampled. 

5. Conclusion and future work 

We have presented the detection problem of subpixel objects embedded in additive Gaus- 
sian noise. Subpixel location and signal amplitude are assumed to be unknown. Unknown 
subpixel location has a great influence on detection performance in the aliased case while 
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Fig. 10. Average mean square errors (MSE) of position estimators in the Gaussian 
white noise case with common sensor design on the left (r c = 2.44) compared to a 
correctly sampled optics on the right (r c = 0.5). 



conventional matched filter neglects it. Thus, we derived four types of improved detectors 
from the likelihood ratio: the GLRT, the ELRT, the ALRT and the SM-GLRT. We have il- 
lustrated their performance in comparison with the more classical GPMF. Numerical results 
for both white and correlated noise cases show that the ELRT, the ALRT and the GLRT are 
competitive whereas the SM-GLRT does not reach the same quality but slightly improves 
the performance of the GPMF too. The ALRT seems to be a good trade-off since it is not 
as computionnally demanding as the ELRT and the GLRT. Moreover the performance gain 
proves to be only moderate in the case of unaliased optics. This conclusion has important 
consequence in sensor design: it suggests that the popular design of a pixel covering exactly 
the main lobe of the Airy disk is not optimum for point object detection. Future work con- 
sists in studying the robustness of these detectors to real data and the way we can take into 
account non Gaussian distributions of background noise. As far as the position estimation 
problem is concerned, we have demonstrated prospective gains that must also be confirmed 
on more realistic data. 
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This paper addresses the issue of detecting point objects in a clutter background 
and estimating their position by image processing. We are interested in the specific 
context where the object signature significantly varies with its random subpixel 
location because of aliasing. Conventional matched filter neglects this phenomenon 
and causes consistent loss of detection performance. Thus, alternative detectors are 
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proposed and numerical results show the improvement brought by approximate and 
generalized likelihood ratio tests in comparison with pixel matched filtering. We 
also study the performance of two types of subpixel position estimators. Finally, 
' we put forward the major influence of sensor design on both estimation and point 

object detection performance. © 2009 Optical Society of America 

1. Introduction 

We tackle the problem of subpixel object detection in image sequences which arises for in- 
stance in infrared search and track (IRST) applications. In this context, the target signature 
is proportional to: 

/•i+0.5 /-i+0.5 

s e [i,j]= / / h a (u - e u v - e 2 ) dudv. (1) 

Ji-0.5 Jj-0.5 

8 e [i,j] represents the percentage of light intensity at pixel e = (ei, e 2 ) refers to the 
object random subpixel position and h a is the optical point spread function (PSF). Ac- 
cording to common sensor design, the energy of the signal component s = as e is almost 
concentrated on a single pixel. However, contrary to the amplitude a which is unknown too, 
dependence on the location parameter e is highly nonlinear. Its influence in our application 
is rather significant because of aliasing and unless a velocity model is available, object sub- 
pixel position is hardly predictable from frame to frame. Actually, common sensor design 
leads to an image spot downsampled by almost a factor 5. We can see on Figured] the en- 
ergy loss at central pixel according to subpixel location and the random change in spatial 
pattern due to aliasing. This phenomenon has a major impact on detection performance as 
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Fig. 1 . Examples of image spots for different cross-marked subpixel positions (win- 
dows of size 5x5 pixels). Sensor design parameter r c is set to its common value of 
2.44 (see section HJ). 



shown thereafter. To our knowledge, this pitfall has not been addressed yet in the litera- 
ture. A prevailing opinion stands that there is no signature information in subpixel objects. 
Indeed, the different authors dealing with small object detection have concentrated on clut- 
ter removing, 13 multi- or hyper- spectral fusion 4 5 and multiframe tracking methods. 6-8 We 
focus here on the processing of a single frame. In section [2l we formulate the detection 
problem in the classical model of a signal in additive Gaussian noise [9, ch.2-4]. When the 
signal is deterministic, Neyman-Pearson strategy yields the conventional matched filter. In 
the present case, the signal from the target depends on unknown parameters and we have 
to deal with a composite hypothesis test. A common procedure is given by the generalized 
likelihood ratio test. But the " nuisance " parameters a and e can also be considered as ran- 
dom variables with known distributions (some a priori density functions in the Bayesian 
terminology), then the straightforward extension of the likelihood ratio test is to integrate 
the conditional distribution over a and e. When modelling the signal component as a sam- 
ple function, we could also think of the class of random signal in noise detection problems, 
which have essentially been studied in the Gaussian case. Unfortunately, considering s e as 
a random vector, its empirical distribution proves to be highly non Gaussian when e is uni- 
formly sampled. For instance, the histogram of the central pixel depicted on Figure|2]shows 
that a Gaussian fit is not satisfactory at all. In section|3l we define more precisely the optical 
system model used in our numerical experiments. We consider both a Gaussian white noise 
and a fractal noise of unknown correlations generated by a standard technique of spectral 
synthesis. Section|4]is devoted to the position estimation problem, i.e. estimation of param- 
eter e. We propose two estimators that take into account the fact that the signal amplitude 
a is also unknown. We demonstrate the performance of these estimators in terms of mean 
square errors. As for the detection problem, we finally illustrate the expected improvement 
in quality brought by a correctly sampled optics compared to common sensor design. 
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Fig. 2. Empirical distribution of the image-spot central pixel s e [0, 0] for a uniformly 
random position e ~ W[-o.5,o.5[ 2 - 



2. Detection problem 

We consider a local detection window sliding across the image. The problem is to decide 
whether an object is present or not at the window central pixel. This is a binary test which 
typically reads as follows: 

H : z = n 
Hi : z — as e + n 

where z is the vector collecting the window data, s = as e is the object response (signal 
vector) and n the additive Gaussian noise. The signature shape is known and deterministic, 
so that s only depends on the two unknown parameters a G R and e e £ = [—0.5, 0.5 [ 2 . 
The noise vector n is supposed to be centered (in practice we first remove the empirical 
mean from the data) with a known or previously estimated covariance matrix R. Thus, if we 
assume that n is independent from s, the following conditional distributions are Gaussian: 

p(z|H ) ~ Af(0,R) 
p(z|Hi,a,e) ~ Af(as e ,R) 

Let first assume that parameters a and e are given. The problem amounts to a simple 
hypothesis test which is to detect a deterministic signal in a Gaussian noise. The Neyman- 
Pearson strategy or likelihood ratio test (LRT) is given by: 



Hi 

p(z|Hi,q,e) > 

p(z|H ) < 
H n 



threshold (4) 



It is equivalent to classical matched filtering which simply compares the statistic aT e (z) = 
a s l e R~ l z with some threshold. 

A. Pixel matched filtering 

The exact object location being unknown in practice, we could assume by default that 
e = e = [0,0], i.e. the object is at the center of the pixel, whereas the true location 



3 



would correspond to e = e*. Thus, the detector which consists in thresholding the pixel 
matched filter (PMF) aT eo (z) is optimum provided e* = e . Otherwise it is mismatched 
and therefore suboptimum. Since conditional distributions of T eo (z) under each assump- 
tion are Gaussian, we easily get the expression of the probalility of detection Pd and of 
false alarm Pf a . Corresponding receiver operating characteristic (ROC) curves for critical 
values of e* are depicted on Figure [3] They clearly show that the PMF performances are 
significantly worse as e differs from e*. But beyond extreme situations (related to a true 
target location between two or four pixels instead of the center), the "mean curve" repre- 
sents the average statistics over uniformly random positions. Compared to the ideal curve, 
we can see that the price paid if one neglects the random location is rather high even at 
favorable signal-to-noise ratio. For a SNR of 15dB and at a Pj a of 1CT 4 , probability of 
detection decreases from nearly 1 to 0.8. 




io~ 1 ° icr s io~" io~* icf s 10° 
Probability of False Alarm 



Fig. 3. Examples of pixel matched filter theoretical ROC curves for e* = e in blue 
(ideal curve), e* ^ e in green and red (worst case), and finally the mean curve in 
black for uniformly sampled e* (SNR = 15dB). 



The object response also depends (linearly this time) on the amplitude a, which is 
generally unknown. Yet, assuming strictly positive amplitude, we see that whatever a > 
0, thresholding aT eQ (z) gives the same ROC curve as thresholding T eo (z). Without any 
assumption on a, a classical solution is to estimate it by maximum likelihood (ML). Indeed 
under the Gaussian noise assumption, the optimum in a for a given e is explicit: 

d(e) = argmax p(z|H 1; a, e) 

= argmin {(2 — as e ) t R~ 1 (z — as £ )} 

siR^z 

and then the " generalized " pixel matched filter (referred to as GPMF) is equal to 

Is* R^zl 2 

a(e )T eo (z) = ' 60 ' . (6) 
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B. Subpixel detectors 

Our aim is to build refined detectors that improve performance of the above GPMF in taking 
into account the variability of the object signature due to its random subpixel location. 
Several solutions may be used. We first recall the most popular one. 

1. Generalized likelihood ratio test 

ML estimation of the two unknown parameters leads to the generalized likelihood ratio test 
(GLRT): 

max (ai£) p(z|Hi,q,e) 
' M) " p(*|Ho) 

(7) 

= > threshold. 

It consists in estimating the amplitude a and the possible object location e by computing: 

e ML = argmax p(z|Hi, a(e), e) 

e<=£ 



arg max 



s{R- 1 z\^ (8) 



then thresholding the estimated filter a ML % ML (z) where d ML = a(e ML ) is given by equation 
©: 

& ML % ML (z) = 1 e - . (9) 



si R-'s, 



2. Exact likelihood ratio test 



In a Bayesian approach, we propose to consider the two unknown parameters a and e as 
realizations of independent random variables with given probability density functions p(a) 
and p(e). Then the optimal procedure is the exact likelihood ratio test (ELRT). 

To compute the density function of data under Hx and to get the likelihood ratio, we 
have to integrate the conditional density p(z|Hi, a, e) over prior distributions of the nui- 
sance random parameters a and e. The likelihood ratio can be expressed as: 

= p(z|Hi) = f e / R p(z|Hi, a, e)p(at)p(e) da de 
[Z) p{z\R ) " p(z\U ) 

Given prior distributions p(a>) and p(e), C(z) is the optimal Neyman-Pearson test whenever 
a and e really satisfy the models p(a) and p(e) . By default we choose a " non-informative " 
prior for a and we adopt a uniform distribution inside the pixel for e, which seems to be 
quite a reasonable assumption for the subpixel target position. So we get: 

C(z) ex / — 1 exp \ i-g- ^ de, (11) 

Unfortunately, because of intricate nonlinear dependence of s e on e, explicit integra- 
tion over e appears to be not tractable and probability distribution of C(z) is not as simple 
as the one of T €Q (z). A quadrature approximation is required to compute C(z) whereas 
derivation of its density requires Monte-Carlo simulations. 
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3. Approximate likelihood ratio test 

In equation (fTTT) . the double integral over e can be approximated up to any desired accuracy 
using some quadrature rule and evaluating the integrand f(e\z) at discrete samples 6 
£ = [—0.5, 0.5 [ 2 . But, for sake of computational efficiency, we propose to use a coarsest 
approximation of the likelihood ratio (ALRT) based on a bidimensional trapezoidal rule 
which only involves the 9 half -pixel positions. 

4. Subspace model 

One alternative to this probabilistic viewpoint can be built on a geometric approach that 
restricts the signal vector s = as e to vary in some P-dimensional subspace, with P lower 
than the vector size. 10 The observed data under H4 are rewritten as: 

p 

z ~ Sa + n = a p s p + n, (12) 
p=i 

where the structural matrix S is formed by P independent vectors s p . Coefficients a p of 
the linear combination are the new parameters that describe the signal variability. Thanks 
to linearity, ML estimation of vector a has an explicit solution (which is identical to the 
least squares estimator): 

a ML = (S'R^S^S'R^z (13) 
and GLRT amounts to threshold the following statistic: 

V{z) = z t Br 1 S(S t Rr 1 S)- 1 S t Rr 1 z. (14) 

Matrix S only depends on e, a being a scale parameter. In practice, it is identified by 
discretizing £ , making a singular value decomposition and retaining the singular vectors 
s p corresponding to the P greatest singular values. We choose P = 1 which gives better 
results than higher orders. Therefore under hypothesis Hi, z ~ aiSi + n and V(z) is 
identical to GPMF with s 1 replacing s eo . 

3. Application to optical imagery 

A. Optical system 

In our application, we can model the imaging system by a diffraction-limited, unaberrated 
optics with circular aperture and incoherent illumination. 11 ' 12 The object signal pattern s e 
is then given by the integration of h a on each pixel (see equation d]), where h is the radial 
point spread function (PSF) defined by the Airy disk: 

2 

, p = Vm 2 + v 2 . (15) 

J\ is the Bessel function of the first kind and r c = v c j v s designates the normalized cut- 
off frequency (y s is the sampling frequency and v c = D/\is the radial cut-off frequency 
defined by the ratio of the lens aperture diameter D over the wavelength A). FigureHJdepicts 
the two-dimensional PSF and a slice along one diameter, as well as their Fourier transform. 



h Q (u,v) 



TV 



Ji{-Kpr c ) 
P 
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Common sensor design uses r c = 2.44 so that the pixel size is equal to the width of the 
main lobe of the PSF However, this implies a downsampling factor of v n /v s = 2r c = 
4.88 (where v n = 2v c is the Nyquist frequency). In the following section, we present 
some numerical results of detection performance considering this classical sensor design. 
Examples of image spots s e have been represented on Figured] for different values of e. 

Remark 1 We have the following property: 




1. 



P 



Fig. 4. Left: radial point spread function h (u,v) on the top and slice along a di- 
ameter on the bottom. Right: corresponding optical transfer function h (u u , u v ) and 
slice along a diameter (r c = 2.44). 



B. Numerical results 

Performances of the five classes of detectors have been compared in terms of ROC curves: 
GPMF, GLRT on a and e, ELRT, ALRT and finally GLRT with the subspace model de- 
noted SM-GLRT. Probabilities of detection and false alarm are deduced from the empirical 
distributions of these statistics under each hypothesis by generating samples of Gaussian 
noise n and uniformly distributed e in £ = [—0.5, 0.5 [ 2 . The amplitude is assumed to be 
unknown but set to a constant value a in the simulations since we have no information 
about a reliable prior distribution p(a). 

We first consider the case of a Gaussian white noise n ~ A/"(0, a 2 ) . The signal-to-noise 
ratio is then defined by: 

fa 2 E\ 

SNR=10 1og 10 (^J (16) 



with E= I (sel*,]}) 2 de 
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aS 0.95 
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SNR = 16.2dB (a/a = 9) 
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SNR= 14.1dB (a/a = 7) 



Fig. 5. Empirical ROC curves in the Gaussian white noise case with common sensor 
design (r c = 2.44). 



For common sensor design (r c = 2.44), the average energy of the image spot is E ~ 0.52. 
The ROC curves are depicted on Figure \5\ for two different SNR. It shows that the GLRT, 
the ELRT (actually, a refined approximation of it) and the coarse approximation ALRT 
give significantly better performance than the SM-GLRT and GPMF. We also see that the 
performance gain is greater for high SNR whereas it tends to be rather small for low SNR 
and low probability of false alarm. Conversely, if the latter detectors are computationally 
cheap, including the ALRT, this is not so for the GLRT and the ELRT, which are much 
more intensive. 

As complementary results, we have tested the five detectors on a fractal background 
image generated by a variant of the ppmforge software 0. The synthesis algorithm depends 
on the auto-similarity parameter H called Hurst parameter and which is set to 0.7 in this 
experiment. The resulting image depicted on Figure [6] is a realistic simulation of a cloud 
scene. The covariance matrix R of this stationary background is estimated by the empir- 
ical correlations on the whole image. We then compute the performance of the different 
detectors for a given target amplitude as illustrated on Figure [7J The ROC curves look quite 

'http://h30097.www3.hp.eom/demos/ossc/man-html/manl/ppmforge.l.html 
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different from the white noise case but we notice again that the GLRT, ELRT and ALRT 
have similar performance and provide a significant detection gain in comparison with the 
GPMF or the SM-GLRT. 



200 
100 



-100 
-200 
-300 



Fig. 6. Simulation of a cloud fractal image of 200 x 200 pixels (Hurst parameter 
H = 0.7). 




10 10 10 10 10 

Probability of False Alarm 



Fig. 7. Empirical ROC curves obtained on the fractal image of Figure [6] for a true 
(but assumed unknown) target amplitude a = 60 gray levels. 



C. Influence of the optics 

Besides the perfecting and evaluation of subpixel detectors, one additional motivation of 
this work is to analyze the influence of aliasing on detection performance. This is the reason 
why we have also tested the detectors on a correctly sampled optics in order to compare 
their performances with those obtained using a common sensor design. In the correctly 
sampled design, the focal plane is sampled at Nyquist frequency (implying a denser sensor 
array or a smaller lens diameter) so that aliasing is suppressed. Parameter r c of the PSF 
is equal to 0.5 and the signal energy is now spread over several pixels (E ~ 0.08). By 
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Fig. 8. Examples of image spots corresponding to a correctly sampled optics (r c = 
0.5) to be compared to those of Figured] 



comparison, Figure [8]presents the examples of image spots corresponding to such a design. 
Detection performances are depicted on Figure [9] on the right for a SNR of 15dB. We see 
that improved detection has just a moderate impact in this situation. The five detectors have 
a quite similar behavior but at the same SNR they perform much better than in the aliased 
case. The gain in Pf a amounts at least to a factor 10 for all the detectors. Such a result 
speaks in favour of using a denser focal plane for point target detection. 




Fig. 9. Empirical ROC curves in the Gaussian white noise case for a same SNR 
= 15dB, with common sensor design on the left (r c = 2.44) compared to a correctly 
sampled optics on the right (r c = 0.5). 



4. Performance of subpixel position estimators 

Up to now we have focused on the detection strategy. In a second step, once a potential tar- 
get is detected on a given pixel, we are also interested in accurate estimation of its subpixel 
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position. Such a problem has already been addressed, in particular for star position estima- 
tion in astronomical applications. 13 Several types of estimators are possible. We consider 
here the maximum likelihood (ML) estimator and following the Bayesian approach intro- 
duced previously the posterior mean (PM). It is important to note that the signal amplitude 
a is also unknown and therefore we have to estimate it or integrate over it. Indeed it is not 
valid to suppose that the amplitude is known in the context of IRST. 

The ML estimator of e is given in equation ([8]) by replacing a with its estimate e. Actu- 
ally, e ML and a ML = d(e ML ) are identical to joint maximum a posteriori (MAP) estimators 
with non-informative priors on the two parameters. 

The PM estimator is defined as: 



J ep(e\E 1 ,z)de (17) 



where the posterior law is deduced from Bayes'rule: 



p(e\H u z) = ri«P';«)*> (18) 
p(*|Hi) 



p( e ) f 

— — / p(;z|Hi, a, e)p(a) da. 
J R 



p(z 

So, we have to integrate over a and then over e. As previously we consider a diffuse a 
priori on R for a and a uniform law on £ for e. We get the following expression in the 
same way as for the likelihood ratio in equation ([111) : 

1 ( I s* 7? _1 zl 2 1 

p(c|Hi, z) oc - exp ' \ — j 1 . (19) 



We have studied the performance of these two estimators in terms of average mean 
square error (MSE). In practice, the optimization or the integration over e are approxi- 
mated numerically by considering a finite discrete grid of 20 x 20 values e k E 8. Given a 
true position e*, bias and variance of an estimator e are estimated thanks to Monte-Carlo 
simulations. We consider the case of a Gaussian white noise and we vary the signal-to- 
noise ratio. Figure [10] on the left compares ML and PM estimators to the pixel estimator 
which assumes by default that the target location is at the center of the pixel (e = (0, 0)) 
and whose MSE is equal to 1/12. At favorable SNR, the two subpixel estimators are far 
better than the default estimator but the gain decreases when the noise becomes important. 
For a SNR of 15dB, the ML yields an error similar to the default estimator while the PM 
notably has a twice smaller error. By comparison, Figure[lO]on the right shows the estima- 
tion performances obtained in the unaliased case (r c = 0.5) for equivalent signal-to-noise 
ratios. ML and PM logically perform better since the signal is correctly sampled. 

5. Conclusion and future work 

We have presented the detection problem of subpixel objects embedded in additive Gaus- 
sian noise. Subpixel location and signal amplitude are assumed to be unknown. Unknown 
subpixel location has a great influence on detection performance in the aliased case while 
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Fig. 10. Average mean square errors (MSE) of position estimators in the Gaussian 
white noise case with common sensor design on the left (r c = 2.44) compared to a 
correctly sampled optics on the right (r c = 0.5). 



conventional matched filter neglects it. Thus, we derived four types of improved detectors 
from the likelihood ratio: the GLRT, the ELRT, the ALRT and the SM-GLRT. We have il- 
lustrated their performance in comparison with the more classical GPMF. Numerical results 
for both white and correlated noise cases show that the ELRT, the ALRT and the GLRT are 
competitive whereas the SM-GLRT does not reach the same quality but slightly improves 
the performance of the GPMF too. The ALRT seems to be a good trade-off since it is not 
as computionnally demanding as the ELRT and the GLRT. Moreover the performance gain 
proves to be only moderate in the case of unaliased optics. This conclusion has important 
consequence in sensor design: it suggests that the popular design of a pixel covering exactly 
the main lobe of the Airy disk is not optimum for point object detection. Future work con- 
sists in studying the robustness of these detectors to real data and the way we can take into 
account non Gaussian distributions of background noise. As far as the position estimation 
problem is concerned, we have demonstrated prospective gains that must also be confirmed 
on more realistic data. 

References 

1. C. D. Wang, "Adaptive spatial/temporal/spectral filters for background clutter suppres- 
sion and target detection," Optical Engineering, vol. 21, pp. 1033-1038, Dec. 1982. 

2. A. Margalit, I. S. Reed, and R. M. Gagliardi, "Adaptive optical target detection using 
correlated images," IEEE Transactions on Aerospace and Electronic Systems, vol. 21, 
pp. 394-405, May 1985. 

3. T. Soni, J. R. Zeidler, and W. H. Ku, "Performance evaluation of 2-D adaptive predic- 
tion filters for detection of small objects in image data," IEEE Transactions on Image 
Processing, vol. 2, pp. 327-340, July 1993. 

4. X. Yu, L. E. Hoff, I. S. Reed, A. M. Chen, and L. B. Stotts, "Automatic target de- 
tection and recognition in multiband imagery : a unified ML detection and estimation 



12 



approach," IEEE Transactions on Image Processing, vol. 6, pp. 143-156, Jan. 1997. 

5. E. A. Ashton, "Detection of subpixel anomalies in multispectral infrared imagery us- 
ing an adaptive Bayesian classifier," IEEE Transactions on Geoscience and Remote 
Sensing, vol. GE-36, pp. 506-517, Mar. 1998. 

6. I. S. Reed, R. M. Gagliardi, and H. M. Shao, "Application of three-dimensional filtering 
to moving target detection," IEEE Transactions on Aerospace and Electronic Systems, 
vol. 19, pp. 898-905, Nov. 1983. 

7. S. D. Blostein and T. S. Huang, "Detecting small, moving objects in image sequences 
using sequential hypothesis testing," IEEE Transactions on Signal Processing, vol. 39, 
pp. 1611-1629, July 1991. 

8. J. M. Mooney, J. Silverman, and C. E. Caefer, "Point target detection in consecu- 
tive frame staring infrared imagery with evolving cloud clutter," Optical Engineering, 
vol. 34, pp. 2772-2784, Sept. 1995. 

9. H. L. Van Trees, Detection, estimation and modulation theory. Parti , New York: Wiley, 
John, 1968. 

10. D. Manolakis and G. Shaw, "Detection algorithms for hyperspectral imaging applica- 
tions," Signal Processing Magazine, vol. 19, pp. 29-43, Jan. 2002. 

11. J. W. Goodman, Introduction a Voptique de Fourier et a I 'holographic Paris: Masson, 
1972. 

12. R. C. Hardie, K. J. Barnard, J. G. Bognar, E. E. Armstrong, and E. A. Watson, "High- 
resolution image reconstruction from a sequence of rotated and translated frames and 
its application to an infrared imaging system," Optical Engineering, vol. 37, pp. 247- 
260, Jan. 1998. 

13. K. A. Winick, "Cramer-Rao lower bounds on the performance of charge-coupled- 
device optical position estimators," Journal of the Optical Society of America A, vol. 3, 
pp. 1809-1815, Nov. 1986. 



13 



